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Abstract. 

We study the purely relaxational dynamics (model A) at criticality in three- 
dimensional disordered Ising systems whose static critical behaviour belongs to the 
randomly diluted Ising universality class. We consider the site-diluted and bond- 
diluted Ising models, and the ±J Ising model along the paramagnetic-ferromagnetic 
transition line. We perform Monte Carlo simulations at the critical point using the 
Metropolis algorithm and study the dynamic behaviour in equilibrium at various values 
of the disorder parameter. The results provide a robust evidence of the existence of 
a unique model-A dynamic universality class which describes the relaxational critical 
dynamics in all considered models. In particular, the analysis of the size-dependence 
of suitably defined autocorrelation times at the critical point provides the estimate 
z = 2.35(2) for the universal dynamic critical exponent. We also study the off- 
equilibrium relaxational dynamics following a quench from T = oo to T = T c . In 
agreement with the field-theory scenario, the analysis of the off-equilibrium dynamic 
critical behavior gives an estimate of z that is perfectly consistent with the equilibrium 
estimate z = 2.35(2). 
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1. Introduction 

Randomly diluted uniaxial antiferromagnets, for instance, Fe p Zni_ p F2 and Mn p Zni_ p F2, 
have been much investigated experimentally and theoretically [TJ [2j [3j H] . For sufficiently 
low impurity concentration 1 — p, these systems undergo a second-order phase transition 
at T c (p) < T c (p = 1). The critical behaviour is approximately independent of the 
impurity concentration and definitely different from the one of the pure system. These 
results have been successfully explained by the field-theoretical (FT) renormalisation 
group (RG), which predicts the presence of a single universality class associated with 
the paramagnetic-ferromagnetic transition that occurs in Ising systems with quenched 
random dilution. Monte Carlo (MC) results have been contradictory for a long time, 
finding model-dependent critical exponents. In [5] this apparent non-universality was 
shown to be an effect of strong scaling corrections. They are slowly decaying due to the 
fact that the leading correction-to-scaling exponent to is quite small: u = 0.29(2) (see 
Appendix A ). 

The analyses significantly gain accuracy when using improved Hamiltonians, for 
which the leading scaling corrections are suppressed for any thermodynamic quantity, 
and improved estimators, which are such that the leading scaling correction is suppressed 
for any model in the same universality class. MC simulations of different improved 
Hamiltonians [HI [7] confirmed that the static critical behaviour is model-independent, in 
agreement with the FT description, and provided accurate estimates of the static critical 
exponents, v = 0.683(2) and n = 0.036(1) [6j[8j[5]. They are in good agreement with the 
FT perturbative results [9] v = 0.678(10) and n = 0.030(3) obtained by the analysis of 
high-order (six-loop) perturbative expansions (similar results are obtained at five loops 
[TO]). The apparent non- universality observed in previous numerical works was mainly 
due to the fact that scaling corrections were neglected. As a consequence, previous 
studies did not really observe the asymptotic critical behaviour and only determined 
effective exponents depending on all parameters of the investigated model. 

In this paper we extend the analysis to the critical dynamics. We consider a purely 
relaxational dynamics without conserved order parameters, also known as model A [llj, 
as appropriate for uniaxial magnetic materials. Experimental results are reported in 
[T2l [131 [Uj. According to the FT RG (see, e.g., [JJ2, [16j E!)' the dynamic behaviour 
should be the same in all RDIs systems, as is the case for the static criticality. 
Moreover, the leading scaling corrections appearing in dynamical quantities should be 
associated with the same RG operators that control the nonasymptotic behaviour of 
static quantities and thus, they should be characterized by the same exponents as in 
the static case, i.e., by uj = 0.29(2) and u>2 = 0.82(8). As a consequence, in the case of 
improved Hamiltonians, leading scaling corrections should also be absent in dynamical 
quantities. Therefore, the most precise estimates of dynamic universal quantities should 
be obtained in improved models, as in the static case. 

Previous MC studies [HI UHl EQl EH E21 E3] of equilibrium and off-equilibrium 
dynamics apparently have not confirmed the FT general predictions. They have mainly 
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focused on the dynamic critical exponent z, which characterizes the divergence of the 
autocorrelation times when approaching the critical point. In most of the cases they have 
found that z is model dependent and have provided estimates which range from z ps 2.1 
to z 2.6, depending apparently on the method, the favoured values of the dilution 
parameter p, whether it is determined from equilibrium or off-equilibrium simulations, 
etc. In [20] [22] the universality of z was verified, obtaining z ~ 2.6, but the leading 
scaling-correction exponent was not consistent with the static one, as predicted by the 
FT approach. Moreover, this result is inconsistent with the perturbative FT estimate 
obtained from analyses of the perturbative expansions [2U [25], [26] [151 123 [28] [29] at two 
and three loops, which suggest z ~ 2.18. 

In this paper we study three disordered Ising systems whose static critical behaviour 
belongs to the 3D RDIs universality class: the randomly site-diluted Ising model 
(RSIM), the randomly bond-diluted Ising model (RBIM), and the ±J Ising model 
along the paramagnetic-ferromagnetic transition line. Their static critical behaviour 
was carefully investigated in [6] [7]. In particular, the value p* of the dilution parameter 
corresponding to an improved model was determined for each of them. We simulate 
these models by using the Metropolis algorithm (with a suitable modification in the 
case of the RSIM and RBIM to avoid ergodicity problems, see Appendix B ), which 
does not satisfy any conservation law, and thus allows us to investigate the model-A 
dynamics. We consider cubic lattices of size L 3 with 8 < L < 64. 

The main purpose is to check whether the dynamic critical behaviour is consistent 
with the FT RG, that is with the existence of a unique model-A universality class for 
RDIs systems. We focus on the dynamic critical exponent z, and determine it in the 
RSIM, the RBIM, and the ±J Ising model. We find that the autocorrelation times 
extracted from the autocorrelation function of the magnetic suspectibility at T c behave 
as 

r = cL z (l + cniT*" + c 12 L~ 2uJ + ■■■ + c 21 L-" 2 + • ■ ■) (1) 

with a universal value of the dynamic exponent z. Moreover, u and u>2 are consistent 
with the static scaling-correction exponents u = 0.29(2) and co> 2 = 0.82(8). We obtain 
the estimates z = 2.355(16), z = 2.335(18), and z = 2.345(17), respectively for the 
RSIM, the RBIM, and the ±J Ising model at p ~ p*. They are in good agreement, 
strongly supporting universality. Results for other values of p, both larger and smaller 
than p*, are consistent with the estimates of z obtained at p ~ p*. We consider 

z = 2.35(2) (2) 

as our best estimate of z for the dynamic model-A universality class of RDIs systems. 
These results confirm the general picture that comes out of the FT analysis. However, 
from a quantitative point of view, our estimate significantly differs from the perturbative 
result z ~ 2.18 at three loops [28] [29]. Apparently, perturbative FT expansions at this 
order are not able to predict accurately the exponent z. 

The exponent z can also be determined by performing off-equilibrium simulations, 
since the approach to equilibrium is controlled by the same FT model [3Ql EU [16]. As 
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a further check of our result (j2J), we have performed off-equilibrium MC simulations of 
the RSIM at p — 0.8, quenching T = oo configurations to T = T c . The results show 
that the relaxation to equilibrium is controlled by the same dynamic exponent obtained 
in equilibrium simulations, i.e. z = 2.35(2). Moreover, the large-time corrections are 
consistent with what is predicted by the FT RG, which relates them to the static leading 
and next-to-leading scaling-correction exponents to = 0.29(2) and u 2 = 0.82(8). Our 
results therefore confirm the FT analysis of the off-equilibrium relaxational dynamics 

[301 ED EE]. 

The paper is organized as follows. In Sec. [2] we define the disordered Ising models 
that are considered in the paper. In Sec. [3] we define the quantities that are measured 
in the MC simulation and discuss the FT predictions. In Sec. H]we report the finite-size 
scaling (FSS) analysis of equilibrium MC simulations of the RSIM, the RBIM, and the 
± J Ising model. In Sec. |5]we study the off-equilibrium relaxational critical behaviour of 
the RSIM, in a quench from T = oo to T c . Finally, we draw our conclusions in Sec. El In 
Appendix A| we refine the estimate of the leading scaling correction exponent, obtaining 



uj = 0.29(2). Some details on the MC algorithm are discussed in Appendix B 



2. Models 



We consider the randomly site-diluted Ising model (RSIM) with Hamiltonian 

H p = - 22 PxPy a x°y ( 3 ) 

<xy> 

where the sum is extended over all nearest-neighbour sites of a simple cubic lattice, 
o~ x are Ising spin variables, and p x are uncorrelated quenched random variables, which 
are equal to 1 with probability p (the spin concentration) and with probability 1 — p 
(the impurity concentration). We also consider the randomly bond-diluted Ising model 
(RBIM) in which the disorder variables are associated with links rather than with sites. 
It is defined by the Hamiltonian 

= ~ X] i*» 0x0 y ( 4 ) 

<xy> 

where the couplings j xy are uncorrelated quenched random variables, which take values 
0,1 with probability distribution 

P(jxy) = pS(jxy - 1) + (1 - P)5(jxy). (5) 

Note that the exchange interaction is ferromagnetic in both models. 

MC simulations [61 [7] have provided strong numerical evidence that the static 
critical behaviour of the RSIM (for 1 > p > p s , where p s is the site-percolation point, 
p s = 0.3116081(13) on a simple cubic lattice [32]) and of the RBIM (for 1 > p > pt,, where 
Pb is the bond-percolation point, pf, = 0.2488126(5) on a simple cubic lattice [33]) belong 
to the same universality class. The most precise estimates of the static critical exponents 
have been obtained by MC simulations: [61 [H [5] v = 0.683(2) and n = 0.036(1). These 
estimates are in good agreement with the perturbative FT results [TO] v = 0.678(10) 
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and 77 = 0.030(3), and with experiments [TJ [2]. Also the leading and next-to-leading 
correction-to-scaling exponents have been computed. Here we shall obtain a precise 
estimate of the leading exponent u, u = 0.29(2), by a combined analysis of the data 
obtained in [6] and those obtained in the present work; see |Appendix A| for details. As 
for the next-to-leading exponent, we quote the FT estimate obtained in [6], LO2 = 0.82(8). 

We also consider the ± J Ising model, defined by Hamiltonian (j3J) with exchange 
interactions j xy which take values ±1 with probability distribution [31] 

P(jxy) = Pttizy - 1) + (1 - p)S{jxy + 1)- (6) 

Unlike the RSIM and the RBIM, the ± J Ising model is frustrated for any p. Nonetheless, 
the paramagnetic-ferromagnetic transition line that occurs in this model for < p < 
1 — pn and p^ < p < 1 also belongs to the RDIs universality class [7]. Here p^ is the 
location of the magnetic-glassy Nishimori multicritical point, which has been recently 
computed in [35J: p N = 0.76820(4). 

In this work we consider a relaxational dynamics without conserved order 
parameters, i.e. the so-called model A. In lattice systems this dynamics is usually realized 
by using the Metropolis algorithm. In the case of the RSIM and of the RBIM however, 
if a sequential updating scheme is used, the Metropolis algorithm with the standard 
acceptance probability P A =min[l, exp(— (3A?i)] is not ergodic and thus it does not 
provide the correct dynamics. An ergodic dynamics is obtained by introducing a simple 
modification which is described in Appendix B In the ±J Ising model we use the 



standard Metropolis algorithm with a sequential updating scheme. In this model the 
specific problem we observed in the RSIM and in the RBIM is not present (note, however, 
that, to our knowledge, a rigorous proof of ergodicity is lacking for this updating scheme; 
this is also the case of the pure Ising model). 

Note that the algorithm with sequential updating does not satisfy detailed balance 
and hence does not strictly correspond to a reversible dynamicsjj] Detailed balance is 
satisfied only if the spins are updated in random order. It is commonly accepted that 
these two dynamics belong to the same universality class: these violations of detailed 
balance are irrelevant in the critical limit. 



3. Autocorrelation times: definitions and critical properties 

We consider the two-point correlation function 

G(x 2 - xi, t 2 - h) = (a(x 1 ,ti) a(x 2 , t 2 )), (7) 

where the overline indicates the quenched average over the disorder probability 
distribution and (• • •) indicates the thermal average. Near the critical point correlations 

I The Metropolis update is obtained from a single-site update. If P z = {pi Z y} is the transition matrix 
for the update of site z, P z satisfies the detailed-balance condition n x Pxy — ir y p y z ] . However, this 
does not imply that the dynamics is reversible. Indeed, if lattice sites are updated sequentially, the 
transition matrix for a full sweep is P sw — p z ^p z ^p z 3 _ _ _ p z " 5 where n is the number of lattice sites. 
P sw does not satisfy the detailed-balance condition since the matrices P Zi , P Zj for nearest neighbours 
Zi and Zj do not commute. For a more detailed discussion, see, e.g., [36] . 
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develop both in space and time. They can be characterized in terms of the equal-time 
second-moment correlation length £ and of an autocorrelation time r. In the infinite- 
volume limit the correlation length £ can be defined as 

i2 _ 1 dG(k,0) 



(8) 

fe 2 =0 



X dk 2 

where G(k,t) is the Fourier transform of G(x,t) with respect to the x variable and 

X = ^G(a;,0) = G(0,0) (9) 

X 

is the static magnetic susceptibility. On a finite lattice with periodic boundary 
conditions, we define £ as 

^ 2 _ G(0,0)-G(g min ,0) 

where g m i n = (2ir/L, 0,0), g = 2 sing/2. To define the autocorrelation time, we consider 
the autocorrelation function A(t) of a long-distance quantity. Then, we define the 
integrated autocorrelation time 

2 t £-^A(0) 2 + ^A(0)' 1 j 

Here t is the Metropolis time and one time unit corresponds to a complete lattice sweep. 

In the critical limit £ and the autocorrelation time r int diverge. If t r = (T — T c )/T c 
and T c is the critical temperature, for \t r \ — > we have in the thermodynamic limit 

£ ~ \tr\~", T int ~ l^r 2 " ~ € Z , (12) 

where z/ is the usual static exponent and z is a dynamic exponent that depends on the 
considered dynamics. 

The correlation function G(k, t) is the quantity of direct experimental interest and 
thus we could take A(t) = G(k,t). However, for the determination of the dynamic 
critical exponent z, it is computationally more convenient to use a different quantity. 
We consider the autocorrelation function of the magnetic susceptibility 

A(t) = (S(0)S(t)) - W, S(t) = ±[j2*(x,t)] 2 . (13) 

X 

Using (iTTj) we could determine the autocorrelation time r int and then, we could use it to 
determine z. However, the determination of this quantity requires the knowledge of the 
large-t behaviour of A(t). Since it is difficult to determine it precisely, r int is unsuitable 
for a high-precision study. We now introduce a new time scale which is particularly 
convenient numerically. Let us define 

n 

T«{t + nl1> S (14) 
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where n is a fixed integer number. A linear interpolation can be used to extend r e ff(t) 
to all real numbers. Then, for any positive x, we define an autocorrelation time t x as 
the solution of the equation 

T x = T CS {XT X ). (15) 

This definition is based on the idea that, if A(t) were a pure exponential, i.e., A(t) = 
Ao exp(— t/r), then r e s(t) = t for all t and thus r x = r for any x. 

Let us now consider the thermodynamic limit with T > T c (high-temperature phase) 
and let us prove that, if the autocorrelation functions decay faster than any power of 
t in the critical limit, then t x behaves as £ z as any "good" autocorrelation time. More 
precisely, we show that T x /T mt is finite and nonzero in the critical limit for any finite x. 
Since A(t) is an autocorrelation function of a long-range quantity, close to the critical 
point it obeys the scaling law 

jjfi = K S ), S = t/r int . (16) 

In the critical limit and for fixed n, we have n/r int — > 0. Thus, we can expand 

r eS (t + n/2) = -r int x j^r [1 + 0{n/r int )] . (17) 

If we now define /Tint, we obtain in the critical limit the equation 

a x = -xf(a x )/ f'(a x ). (18) 

It is a simple matter to show that, if f(S) decays faster than any power of S (S q f(S) — > 
for S — > oo and any q), there is always (at least) one strictly positive solution a x of (TT8"]) I§1 
Thus, we have proved that, for any x > 0, the ratio r x /r int is finite and strictly positive 
in the critical limit. It follows that t x diverges as £ 2 in the critical limit. 

The condition that f(S) decays faster than any power of S is obviously satisfied 
if f(S) decays exponentially, i.e. if f(S) ~ AS a exp(— bS) for large S, where a is 
some exponent. While an exponential decay of the correlations is obvious in pure 
ferromagnetic systems for temperatures T > T c , in the case of random systems some 
discussion is needed. Indeed, in dilute systems one expects a non-exponential relaxation 
for large values of t [37] , due to the presence of rare compact clusters without vacancies 
that are fully magnetized at temperatures that are below the critical temperature of the 
pure system (the same clusters are responsible for the weak Griffiths singularities in the 

§ Proof. The function f(y) is expected to be positive and strictly decreasing, so that f(y) > and 
f'(y) < for any y. Since y q f{y) — > for y — > oo and f(y) > 0, y q f(y) decreases for large values 
of y. Therefore, we have {y q f(y)Y = qy q ^f{y) + y q f'{y) < 0. This implies yf'(y)/f(y) < -q for y 
large enough. Since q can be arbitrarily large, this implies yf'{y)/f{y) — > — oo for y — > oo. To end the 
proof, define h(y) = y + xf(y)/f'(y). For y = 0we have h(0) = xf(0)/f'(0) < 0. For y — > oo, we have 
h{y) = y[l + xf{y)/{yf'{y))] w y — > +oo (here we use the result yf'(y)/f(y) — > — oo for y — > oo). The 
function h(y) is therefore negative for small y and positive for large y. Since it is continuous, h(y) must 
vanish at a finite nonvanishing value of y. 
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high-temperature free energy [38J). For instance, in Ising systems the infinite- volume 
spin-spin autocorrelation function G(x = 0, t) is expected to decay as [37J EHl HOI HI] 

G{x = 0, t) » B exp[-C(ln tf' 2 } (19) 

for t — > oo. In the infinite-volume limit also A(t) may show a non-exponential behavior 
for large £ in the high-temperature phase. However, note that this does not necessarily 
imply that the scaling function defined in (fT6l) decays non-exponentially. On the 
contrary, one can argue [37J that the Griffiths tail ffl9l) becomes irrelevant in the critical 
limit. This is essentially due to the fact that B and C that appear in ( fl9|) are expected 
to be smooth functions of the temperature that approach finite constants as T — > T c . 
Thus, in the critical limit, t — > oo, T — > T c at fixed S, the non-analytic contribution 
simply vanishes!]]] 

In the above-presented discussion, r x represents an infinite-volume autocorrelation 
time determined in the high-temperature phase. A similar discussion applies if we 
consider the FSS behavior. For instance, at T c we have 

§f r f FS s(S), 3 S tL-, (20) 

where L is the lattice size. The function A(t, L) decays exponentially for any L (this 
is rigorously true for an aperiodic dynamics in a discrete spin system). This fact does 
not necessarily imply that f'Fss{S) decays exponentially (a non-exponential behavior 
could occur if the exponential decay sets in for t > t* ~ L z+e , e > 0), though the 
discussion presented above makes this possibility quite unlikely. In any case, if Jfss{S) 
decays faster than any power of S, the previous proof indicates that t x /L z is finite in 
the critical limit for any finite x, and thus t x is a good autocorrelation time. 

Beside the integrated autocorrelation time one can also define an exponential 
autocorrelation time: 

\t\ 



T, 



ex P = - lim -r-r • (21) 
\t\-*co In A(t) 

This quantity is well defined in a finite volume since A(t, L) decays exponentially, but, as 
a consequence of f[T9~j) . it diverges in the infinite- volume limit for all T c < T < T c (p = 1). 
As a consequence, in the infinite-volume limit at fixed temperature, r x diverges as 

j This phenomenon can be easily understood if one imagines A(t) to have the form 
Ait; a exp(— B\t^ z ) + A2 exp(— i?2(ln t) 3 / 2 ). The first term is the critical contribution, while the second 
one is the non-exponential Griffiths tail. The second term dominates for t ^S> t* , where t* is the value 
of t at which the two terms have the same magnitude. In the critical limit we have t* ~ £ z (ln£) 3 / 2 . 
Since the critical limit is taken at t/£ z fixed, the relevant quantity is t* , which diverges as (ln£) 3 / 2 , 
as T — > T c . This means that, for any fixed value of S = t/T mt ~ t/!; z , sufficiently close to the critical 
temperature, t always satisfies the condition t <^ t* , i.e. belongs to the region in which the non- 
exponential tail (I19|) is negligible. These considerations also indicate that one should limit oneself to 
times t <C t* in studies of the infinite- volume critical behavior in the high-temperature phase. Therefore, 
one should always choose x so that t x <C t* for all considered systems. Otherwise, the extrapolated 
critical behavior would be incorrect. 
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x — > oo. However, the decoupling of the non-exponential tail in the critical limit implies 
that 

lim lim ^ (22) 

is finite and related to the decay rate of f(S) for large S. Of course, the two limits in 
fl22|) cannot be interchanged. 

On a finite lattice of size L 5 7"e X p is always well defined. Nonetheless, this does 
not imply that r exp is a good autocorrelation time. On the contrary, at T = T c we 
expect T e ^ p /L z to diverge as L — > oo. Indeed, for each L, r exp (L) is always given by 
the decay rate of the autocorrelation function for the slowest sample, however small is 
the amplitude of this contribution to the autocorrelation function (note that, for finite 
values of L, the disorder average is a finite sum). As a consequence, r exp (L) is the 
exponential autocorrelation time for a pure Ising system in the low-temperature phase, 
which is expected to increase faster than any power of L, as L — > oo [if tunnelling events 
dominate r exp (L) ~ exp(2aL 2 )]. Therefore, r cxp /L z — > oo as L — > oo. The irrelevance 
of the Griffiths phenomenon in the critical limit should however imply that 

lim i im T *Pc,L) (23) 

X— HX L^OO _L Z 

is finite and related to the decay rate of fFss(S). This is the finite- volume analogue of 
§2}. 

In the definition (TH1) the integer n can be taken arbitrarily. However, the 
asymptotic critical behaviour is observed only if n <C r- mt , see (ITTj) . Therefore, in 
practice n should not be too large. It is also convenient to take n not too small, since 
this avoids computing A(t) too frequently in the MC simulations. Note also that, when 
n decreases, the errors on r e ff(t) increase since A(t) and A(t + n) are close. The effect is 
however small, because of statistical correlations that also increase as n decreases. In our 
work we have always considered values of n much smaller than r x (typically n < t x /20) 
and we have verified that the estimate of the autocorrelation times are independent of 
the chosen (small) value of n. 

Definition (|14p provides an effective exponential autocorrelation time at a finite 
time scale. In the same spirit, one can also define truncated integrated autocorrelation 
times. Define 

for any integer k, and I(t) for any real t by linear interpolation. Then, we can define 
an autocorrelation time r x ^ nt as the solution of the equation 

T x,int = I( XT x,mt)- (25) 

For any x, this definition provides a good autocorrelation time, which converges to r int 
for x — > oo. This definition is similar to that proposed in [32]; note, however, the 
completely different spirit in the two definitions. In [32] the method was proposed as a 
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practical self-consistent method for the determination of Ti nt and for this reason x had 
to be large (in practice x was usually taken between 5 and 10). Instead, if one is not 
interested in determining T- mt but only in computing z, x can be taken at will. 

In this paper we compute the exponent z from the volume dependence of an 
autocorrelation time at the critical temperature. Including scaling corrections, we expect 
a behaviour of the form 

r = cL z (1 + c u L-" + c 12 L~ 2uJ + ■■■ + c 2l L~^ + •••), (26) 

where u and u 2 are the leading and next-to-leading critical exponents. As in [U [7] 
we also consider the dynamical behaviour at a fixed value of a renormalized coupling 
constant. Also in this case autocorrelation times behave as in (|26l) . 

In order to determine z it is crucial to have some knowledge of the correction-to- 
scaling exponents that appear in (126]) . RG predicts that the static correction-to-scaling 
exponents also occur in dynamic quantities. For instance, if x = G(0, 0) behaves as 
aL 2_,? (l + enL~ Wstat ) at criticality for L — > oo, then a correction term decaying as 
£-u>stat j g a j go eX p ec ted in G(0,t) for any t ^ 0. However, dynamics gives also rise to 
new scaling corrections and they may decay slower than the static ones (for instance, this 
occurs in the model-C dynamics, see Sec. [6]). In this paper we make the assumption that 
no new scaling corrections with exponent less than u 2 = 0.82(8) appear, as indicated 
by the FT description of the model-A dynamics. As we shall see, this will be confirmed 
by our numerical analysis. Thus, in (1261) u and u 2 should be identified with the static 
scaling-correction exponents. 

In our analysis, we make use of improved models, which are such that the leading 
scaling correction with exponent u vanishes. Since ratios of leading scaling-correction 
amplitudes are universal (both in static and in dynamic correlation functions), this 
cancellation also occurs in dynamic quantities. Improved models have been determined 
in p[7j: the RSIM at p* = 0.800(5), the RBIM at p* = 0.54(2), and the ±J Ising model 
at p* = 0.883(3) are improved. In these models the scaling corrections proportional 
to L~ kw vanish, so that the leading correction-to-scaling exponent is u 2 . Therefore, 
numerical studies of improved models are expected to provide the most precise estimates 
of universal quantities. Of course, this is true only if the usual model-A FT description is 
correct; otherwise, there could be corrections with a new dynamic exponent Ud yn < u 2 , 
which do not cancel and may give rise to large corrections even in models that are 
improved for static quantities. A stringent check of this picture should be the fact that 
the three different improved models we consider give consistent results. 

4. Equilibrium estimate of the dynamic critical exponent z 

4-1. Monte Carlo simulations 

We perform MC simulations of the RSIM, the RBIM, and the ±J Ising models for 
various values of p, close to the critical temperature on cubic lattices of size L 3 with 
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Table 1. MC estimates of t x (L) for the RSIM at p — 0.8 and for various values 
of x at (3 = 0.285744. For x = 1 we also report estimates of t x {L) extrapolated to 
f3 c = 0.2857431(3). We also report the value of n that enters in the definition (fl4|) . 



L 


n 


x = 0.6 


x = 1 


x — X citj l3 c 


x = 1.5 


x = 2 


8 


1 


7.311(5) 


7.946(8) 


7.946(8) 


8.342(15) 


8.535(25) 


12 


2 


18.016(10) 


19.827(17) 


19.826(17) 


20.88(3) 


21.35(5) 


16 


2 


34.783(20) 


38.42(3) 


38.42(3) 


40.57(6) 


41.57(10) 


24 


4 


88.47(5) 


98.21(8) 


98.20(8) 


103.78(14) 


106.51(24) 


32 


6 


172.25(9) 


191.64(16) 


191.61(16) 


202.8(3) 


207.9(5) 


48 


16 


442.4(3) 


494.0(6) 


493.8(6) 


523.1(1.0) 


538.5(1.7) 


64 


30 


864.0(1.1) 


966.4(2.0) 


966.0(2.0) 


1024(3) 


1052(7) 



Table 2. MC estimates of t x=1 (L) for the RSIM at p = 0.85 and (3 C = 0.2661561(5), 
p = 0.8 and f3 c = 0.2857431(3), and p = 0.65 and C = 0.370168(2), and for the RBIM 
at p = 0.7 and & = 0.326710(3), and p = 0.55 and C = 0.432291(2). 



L RSIM p = 0.85 RSIM p = 0.8 RSIM p = 0.65 RBIM p = 0.7 RBIM p = 0.55 



8 


7.595(7) 


7.946(9) 


10.343(10) 


9.410(19) 


12.853(14) 


12 


18.322(13) 


19.826(17) 


30.79(3) 


22.746(22) 


33.30(3) 


16 


34.564(24) 


38.42(3) 


67.55(6) 


42.64(4) 


65.41(5) 


24 


84.94(5) 


98.20(8) 


204.66(22) 


103.45(6) 


169.03(11) 


32 


161.15(8) 


191.61(16) 


447.7(7) 


193.56(10) 


331.38(21) 


48 


398.0(4) 


493.8(6) 


1326(3) 


468.5(5) 


853.2(1.0) 


64 


756.6(1.4) 


966.0(2.0) 


2846(12) 


874.5(1.9) 


1676(5) 



L < 64 and periodic boundary conditions. We use the Metropolis algorithm with 



multispin coding as described in |Appendix B 



For each lattice size we consider N s disorder samples, with iV s decreasing with 
increasing L, from N s m 64 x 10 5 for L = 8 to N s w 64 x 10 4 for the largest lattice L = 64. 
Note that these numbers of samples are much larger than those typically considered in 
previous numerical studies. For each disorder sample, we thermalize the system by 
using a mixture of Metropolis and Wolff cluster updates in the case of the RSIM and 
of the RBIM, while in the case of the ±J Ising model we only used the Metropolis 
algorithmic Then, at equilibrium, we perform runs of approximately 20r Metropolis 

The presence of rare disorder instances characterized by large compact clusters with no vacancies — 
those that give rise to the Griffiths tail — might be a serious problem for the thermalization if only the 
Metropolis algorithm is used. If a fixed thermalization schedule (independent of the disorder sample) is 
employed, the system may be thermalized on average, but in a few rare cases the sampling may begin 
much before the equilibrium state has been reached. However, the considerations presented in Sec. [3] 
indicate that these contributions are irrelevant for the critical behavior. Moreover, their probability is 
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Table 3. MC estimates of t x —\{L) for the ±J Ising model at various values of p and 
at fixed Z/L = 0.5943. 



L 


±J lsp = 0.83 


±J Is p = 0.8 


!83 ±J lsp — 0.9 


8 


10.260(22) 


6.507(3) 


6.035(8) 


10 




10.918(5) 




12 


29.73(5) 


16.590(7) 


14.882(10) 


14 




23.705(9) 




16 


63.94(11) 


32.300(17) 


28.471(16) 


18 




42.353(23) 




20 




54.13(3) 




22 




67.47(4) 




24 


188.7(3) 


82.53(8) 


70.84(6) 


28 




118.02(17) 




32 


407.1(8) 


161.50(20) 


135.37(19) 


40 




271.1(4) 




48 




413.5(7) 


336.9(6) 


56 




592.9(1.6) 




64 


2540(15) 


813(3) 


652(3) 



sweeps, where r is the typical autocorrelation time. The averages over disorder are 
affected by a bias due to the finite number of measures at fixed disorder [331 [6] . A bias 
correction is required whenever one considers the disorder average of combinations of 
thermal averages. We use simple generalizations of the formulas reported in App. B of 
[6] Q Errors are computed from the sample-to-sample fluctuations and are determined 
by using the jackknife method. 

We considered the RSIM at p — 0.8, 0.65 (which are the same values considered in 
[6]) and also at p — 0.85. For the RBIM we worked at p — 0.7, 0.55. These runs provided 
new data for the static quantities that were merged with the old ones [6] and with the 
results obtained in some additional cluster MC simulations at the largest lattices. They 

quite low. For instance, in the RSIM at p — 0.8, for L = 64 the probability of a cube of size I = 5 
(I — 6) without vacancies is of order 10~ 7 (10~ 16 ), which should be compared with 1/N S w 1.6 • 10~ 6 . 
+ In App. B of [5] we discuss the case of uncorrelated data. In our case correlations are relevant and 
thus we must somehow modify those expressions. For instance, in order to compute (B) 2 , we use 

^ N a N m /2-k N m 

l\ s {l\ m Ik) a=1 , =1 J=Nm/2 +k 

where B is a generic observable, B a ^ are the corresponding MC estimates, N s is the number of samples, 
N m is the number of measures in equilibrium for each sample, and k a suitable number. We have usually 
taken k w 2t. The bias is of the order r 2 C , Bs(2fc)/Ar^ l , where r is the integrated autocorrelation time 
of the variable B, and Cbs(^) the corresponding autocorrelation function. Similar expressions are used 
in other cases. 
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Figure 1. The effective exponent z e g(L) vs L~ u with u = 0.29 for the RSIM at 
p = 0.85,0.8,0.65, as obtained from r x= i(£). The dotted lines correspond to the 
estimate z = 2.355(16) obtained by using the RSIM data at p = 0.8. 



allowed us to obtain a new estimate of ui (see Appendix A ) and new estimates of (3 C . 
Repeating the analysis presented in [6] we obtain (3 C = 0.2857431(3), (3 C = 0.370168(2), 
and p e = 0.2661561(5) for the RSIM at p = 0.8, 0.65, 0.85, respectively, (3 C = 0.432291(2) 
and p c = 0.326710(3) for the RBIM at p — 0.55, 0.7, respectively. For each p and P we 
usually considered two values of P very close to P c and determined the autocorrelation 
times at p c by linear interpolation. For the RSIM at p — 0.8, runs were performed 
at P — 0.2857440 and subsequently extrapolated at p c = 0.2857431 (see below). For 
the ±J model we did not perform additional simulations and used the results of [7j. 
They allowed us to determine t x (L) at £/L = 0.5943. In the case of the RSIM we also 
determined t x (L) at fixed £/L = 0.5943. The results are very similar to those obtained 
at T c , and therefore we do not consider them in the following. 

Estimates of t x (L) for the RSIM at p — 0.8 are reported in Table [TJ In the table 
we report the data at P = 0.2857440 and, for x = 1, also the extrapolations at p c . 
Note that the correction due to the small change in P is significantly smaller than the 
statistical error. Estimates of t x (L) for x = 1 at T c for all models and several values of 
p are reported in Tables [2] and [3j 



4.2. Results for the RSIM 

In order to determine z, we define an effective exponent 

z eS (L) = — , 

which, for T = T C and L — >• oo, behaves as 

z eS (L) = z + e X iL~" + e 12 L~ 2LV + ■■■ + e 21 L- W2 + 



(27) 



(28) 
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Figure 2. The effective exponent z e s(L) and the corresponding improved quantities 
z im (L) and z im2 (L) defined in ([52]) versus L~" 2 with lu 2 = 0.82. Results for the RSIM 
at p — 0.8. The effective exponents are obtained by using t x (L) with x = 1. The 
dotted lines correspond to the results of fits to z + aL~ - 82 with L min = 12. 



see ( 1261) . In Fig. Q] we show z eS (L) as obtained from the estimates of t x= i(L) for the 
RSIM at p = 0.85, 0.8, 0.65, reported in Table [2J The raw data show significant scaling 
corrections and it is far from clear that their limit for L — > oo is independent of p. 

In the following we present a detailed analysis of the MC data for the RSIM. First, 
we analyse the data at p — 0.8. If the FT description is correct, we should observe a 
fast convergence to the infinite- volume limit, with corrections proportional to L - ^ 2 . The 
results presented in Sec . 14 . 2 . H confirm this prediction. In particular, there is no evidence 
of a correction-to-scaling exponent smaller than uj 2 in dynamic quantities. These results 
support the general FT scenario which predicts that the two leading correction-to- 
scaling exponents are the static ones us and u)%. Then, we assume the FT scenario and 
perform a consistency check, verifying that the large differences observed in Fig. [T] can 
be explained by scaling corrections. To make the check more quantitative, we introduce 
an improved estimator for the exponent z (we use here the same strategy employed in [6] 
for static quantities) and show that it converges to the same value obtained for p = 0.8 
with the expected scaling corrections. This allows us to confirm universality, i.e. the p 
independence of the dynamic critical behaviour. 

4-2.1. Analysis at p = 0.8. Let us first analyse t x (L) for the RSIM at p = 0.8. If 
the standard FT description of the model-A dynamics holds, the static correction-to- 
scaling exponents are the most relevant ones. Since the RSIM at p — 0.8 is improved, 
the 0(L~ kuj ) scaling corrections are suppressed and therefore we expect the dominant 
scaling corrections to be proportional to L~ w ' 2 with uj 2 = 0.82(8). In Fig. [2] we plot 
z e g(L) as obtained from t x= i(L) versus L -0 82 . The data with L > 12 clearly fall on a 



Relaxational dynamics in 3D randomly diluted Ising models 



15 




Figure 3. Estimates of the dynamic exponent z for the RSIM at p — 0.8 obtained 
from fits of t x (L) to ([50)1 and of the effective exponents to z + aL~ e . We report results 
corresponding to different values of L m i n : L m - m — 8, 12, 16, 24 (some data are slightly 
shifted along the x axis to make them visible). The dotted lines correspond to the 
estimate ([55]), z = 2.355(16). 



line. To determine z we assume t x (L) to behave as 

t x (L) = cL z (1 + c 2 L- e ) (29) 
for L —* oo, and perform fits of the form 

\nr x (L) = zlnL + lnc + c 2 e- £lnL , (30) 

with e = 0.82, 0.74, 0.90, which correspond to u 2 = 0.82(8). Results for x = 0.6, 1, 1.5, 2 
are shown in Fig. [3] versus L m ; n , the smallest lattice size used in the fit. They are 
independent of L min for L min > 12, with x 2 /DOF < 1 (DOF is the number of degrees 
of freedom of the fit). For example, for x = 1 and e = 0.82, we obtain z = 2.357(4) 
and c = 0.0525(10) for L min = 12, and z = 2.356(6) and c = 0.0526(15) for L min = 16. 
For x = 0.6,1.5,2, e = 0.82, and L min = 12, we obtain z = 2.354(3), z = 2.356(7) 
and z = 2.358(13). One can also estimate z by fitting z e g(L) to z + e2iL~ UJ2 . If we 
determine z e g(L) from r x= i(L), we obtain z = 2.357(4) for L min = 12 and z = 2.357(6) 
for L min = 16. All results are perfectly consistent. From these analyses we obtain the 
estimate 

z = 2.356(6) [3], (31) 

which is the result of the fit of t x= i(L) with e = 0.82 and L m \ n = 16. The error in 
brackets gives the variation of the estimate as u>2 varies within one error bar. 

In the above-reported determination we have implicitly assumed that the RSIM at 
p = 0.8 is exactly improved so that there are no leading scaling corrections. However, 
p* is only known approximately and thus some residual 0{L~~ W ) scaling correction are 
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Figure 4. Estimates of z e g(L), as obtained from t x (L) for several values of x, for the 
RSIM at p = 0.8. The dotted lines correspond to the result z = 2.355(16), see ([55]) . 



still present. To determine their relevance, we exploit the fact that ratios of leading 
scaling-correction amplitudes are universal and use the bound [6] 

\co,n{p = 0.8)/co,ii(p = 0.65)| < 1/30, (32) 

which holds for any quantity O, be it static or dynamic, computed in the RSIM at 
p = 0.8 and p = 0.65 (cqii is the amplitude of the L~ w correction appearing in the 
large-L behaviour of O). Bound (|32|) shows that r x (p = 0.8; L) 1+k r x (p = 0.65;L)~ fc is 
exactly improved (the leading correction proportional to L~ w exactly cancels) for some 
k satisfying \k\ < 1/30. Thus, an upper bound on the systematic error due to the 
scaling corrections is obtained by analyzing 

r x (p = 0.8; L) 1 ^ 30 x r x (p = 0.65; L)^ 30 , (33) 

instead of r x (p = 0.8; L). The estimate of z varies by ±0.008, which represents the 
systematic error due to the residual L~ u corrections. The final result is therefore 

z = 2.356(6)[3]{8}. (34) 

The above-presented analysis shows that the estimates of z obtained by using t x 
with different values of x are perfectly consistent, as of course should be expectedo 
There is therefore little advantage in considering many values of x and it is simpler to 
restrict the analyses to a single x. We wish to choose it in such a way to minimize scaling 
corrections and statistical errors. As is clear from Tabled] statistical errors decrease with 
decreasing x. In Fig. H]we show z e s(L) as computed from r x for different values of x. 
Scaling corrections decrease with increasing x and are essentially independent of x for 
x ^ 1. Thus, a good compromise between small statistical errors and small scaling 
corrections is obtained by taking x neither too small nor too large. We have thus chosen 

* The consistency of the estimates shows also that the potential problems due to the Griffiths tail do 
not occur at the values of x and L we consider here. 
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x — 1. The quantities that are analysed in the following sections are always obtained 
from t x=1 (L). 

4-2.2. Analysis for p = 0.65 and p = 0.85. Let us now consider the RSIM at p — 0.65 
and p = 0.85. Since the model is not improved we must include corrections with 
exponent uj and 2uo at least, i.e. consider correction-to-scaling terms proportional to 
L~ w Pd L~ ' 29 and L~ 2uJ « L~ a58 , which decrease slower than the leading correction 
term L~ W2 = L~ 82 occurring in improved models. Assuming this type of corrections, 
we fitted t x= i(L) with the ansatze 

t x=1 {L) = cL z , (35) 



t x=1 (L) =cL z (l + c u L-) (36) 



and 



t x=1 (L) = cL z (1 + c xl L~" + c 12 L~^) , (37) 
fixing uj = 0.29. 

Let us first discuss the case p = 0.65. Fits to fl35|) give x 2 /DOF « 1 (DOF is 
the number of degrees of freedom of the fit) starting from L min = 32. For L min = 32 
we obtain z = 2.671(5). Fits to (|36|) give x 2 /DOF « 1 starting from L min = 24. For 
-^mm = 24 we obtain 2 = 2.46(2). Fits to (|37l) have x 2 /DOF close to one already for 
£mm = 12. For L min = 12 and L min = 16 we obtain z = 2.31(3) and z = 2.23(6), 
respectively. 

The same analysis can be repeated for p = 0.85. If we consider the smallest L min 
corresponding to x 2 /DOF close to 1 for each fit ansatz, we obtain z = 2.229(1) (fit to 
(I35D, L min = 24), z = 2.33(3) (fit to ([36]), L min = 12), and z = 2.12(2) (fit to (ETJ), 
-^min = 8). Again, the results of fits to (j3Tjl vary significantly with L min : for L min = 12 
we obtain z = 2.21(8). 

Using the simple power-law ansatz ( 1351) one obtains results that apparently indicate 
non- universal, p-dependent values of z. Including the expected corrections to scaling 
the results for the critical dynamic exponent z change rather dramatically, indicating 
that scaling corrections play a crucial role in the analysis. However, since the results 
obtained for z depend strongly on the number of correction terms included in the fit and 
also on the minimal lattice size L m i n , we cannot obtain a direct accurate estimate of z 
at these values of p. Analogously, it is not possible to include the additional correction 
term C2iL~ UJ2 ~ L~ 82 , which was important for the analysis at p — 0.8 (in this case we 
should also consider the equivalent correction L~ 3uj ~ L~ a87 ). For these reasons, we do 
not quote a final result for z at p = 0.65 and 0.85. 



4-2.3. Correction-to-scaling amplitudes. Here we assume that the value of z is 
universal, i.e. that it does not depend on p. Based on this assumption, we compute 
amplitude ratios that involve the correction amplitude en defined in (1281) and verify 
that these ratios do not depend on the chosen value of p. This provides a consistency 
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Figure 5. Difference Az(p; L) = z c g(p; L) — z e g(p = 0.8; L) versus L ^ with u> = 0.29 
for the RS1M at p — 0.85, 0.65. The dotted lines are the results of fits to aL~ u in the 
case of p = 0.65, and to aL^ + bL- 2uJ for p = 0.85. 



check that the dynamic universality class is independent of p. This type of analysis 
is equivalent in spirit to an analysis in which data at different values of p are fitted 
together assuming the same dynamic exponent z. For instance, this is what was done 
in [5]. There are, however, two significant differences: first, we use the static correction- 
to-scaling exponents (this allows us to consider the leading and the subleading scaling 
correction); second, we verify that the amplitudes of the leading scaling correction satisfy 
the constraints imposed by the RG, i.e. we verify the universality of the amplitude ratios. 
For this purpose we consider 

Az{p; L) = z cS {p; L) - z eS (p = 0.8; L). (38) 

For L —>■ oo it behaves as 

z cS {p; L) - z cS {p = 0.8; L) w e 11 L^ + e 12 L~ 2 " + ■■■ + e 21 L~^ + ■■■ (39) 

if the dynamic critical behaviour does not depend on p. Since the RSIM at p — 0.8 is 
approximately improved, we have en ~ for p = 0.8, so that 

en»eu(p). (40) 

In Fig. [5] we show the difference (|39l) as obtained from the available data. Fits of 
Az(p; L) to aL~ u and aL~ u + bL~ £ with e = 2u, u>2 provide estimates of en. We obtain 
en = 0.9(2) at p — 0.65 and en = —0.55(15) at p — 0.85. As expected, corrections have 
opposite sign in the two cases and are quite significant at the present values of L. Note 
that at p = 0.85 only fits with two corrections give a reasonable x 2 > indicating that at 
least two correction terms must be taken into account. 
Then, we consider the static quartic cumulants 
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Figure 6. Ratios S^fe-k) and Sd{p;L) defined in (US")) for p = 0.85,0.65, versus 
L~ u , w = 0.29. 



where Hk = ( (XL <7x ) fe )> a ^ fi xe d £/L = 0.5943 — we call them U22, U4, and Ud, 
respectively. For L — > 00 they behave as 

U# = U* # + c #<n L-" + ..., (42) 

where [6] U% 2 = 0.148(1), Ul = 1.648(3), and U* d = 1.500(1). The ratios of the leading 
scaling-correction amplitudes are universal. In the case of U22 and Ud, we have [B] 

C22,ll 



-0.44(3). 



(43) 



(44) 



Cd,n 

Analogously, the ratio 

= en 

is expected to be universal if the dynamic universality class is independent of p. The 
ratios (|44p can be directly estimated by considering 

z eS {p; L) - z cS (p = 0.8; L) 



S#(p; L) 



2 u s # +b 1 L~ u '+b2L 



U # (p;2L)-U#(p = 0.8;2L) 

In Fig. [6]we show S22 and Sd for p = 0.85, 0.65. At p — 0.65 a fit of the data with L > 16 
to a + bL^ (uj = 0.29) gives s 2 2 = 9(1) and s d = —4.5(5). At p — 0.85, the same fit gives 
S22 = 7.4(9) and Sd = —4.2(5). The agreement is satisfactory, taking also into account 
that the errors do not take into account several sources of systematic uncertainty. The 
approximate p-independence of the ratios s 2 2 and Sd represents a nontrivial check that 
the dynamic universality class is independent of p. Assuming universality, we obtain for 
the RSIM 



(45) 



s d = -4.5(5), s 22 = 9(2). 
Note that these ratios are consistent with Sd/s22 = s c 



(46) 



-0.44(3), cf. (PD. 
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It is interesting to note that the scaling corrections occurring in t x (L) are 
significantly larger than those occurring in static quantities. For instance, we have 

en = s#Upn2 ^ f -5(1) for U 22 
(c#, n /U* # ) ~ 2--1 ~ \26(3) for U d , ' " 

where Cn is defined in (J26J). 

4-2.4- Improved estimators and universality. In the estimate of z obtained at p — 0.8 
in Sec. 14.2.11 cf. (!34l . the residual 0(L _aJ ) scaling corrections are an important source of 
error. These corrections can be significantly reduced by considering improved estimators 
[6]. The estimate of the universal ratio Sd obtained in Sec. 14.2.31 allows us to define 
improved quantities with smaller L _UJ scaling corrections. Let us consider the quantities 

Z 1 (r;L)=z eS {L) (U d {2L)/Uff , 

Z 2 (q; L) = z cS (L) + q(U d (2L) - U* d ). (48) 

For L —>■ oo they behave as 

Zi(L) = z + f n L-" + f 12 L~ 2 " + ■■■ + f 2X L-"> + ■■■. (49) 

where the correction-to-scaling amplitudes depend on r or q. Then, we determine r* 
and q* such that /n(r*) = fn(q*) = 0. An easy calculation gives 



VU* A s d 



q* = -2 u s d . (50) 



z 

Note that r* and q* are expressed in terms of universal quantities and thus Zi(r*;L) 
and Z 2 (q*; L) are improved in any model in the same dynamic universality class. Using 
Sd = -4.5(5), z = 2.36(2), U* = 1.500(1), we obtain 

r* = 3.5(5), q* = 5.5(7). (51) 

The error is mostly due to the error on s d . For r = r* and q = q* the scaling corrections 
are proportional to L~ 2u) with 2uo = 0.58(4). In the following, we define improved 
estimators by taking the central values of the estimates (|5T|) : 

z im (L) ee Z x {r = 3.5; L), z im2 (L) = Z 2 (q = 5.5; L). (52) 

One may define analogous improved operators by using U 22 instead of U d . Those 
defined in terms of U d are more convenient because U d is known with better numerical 
precision. Since r* and q* are known only approximately, z- im {L) and Zi m2 (L) still have 
L~ w corrections. Taking into account the uncertainty on the estimates of r* and q*, we 
obtain the bound 

|/u/en| £ 1/6, (53) 

i.e. the leading scaling correction in Zi m (L) and z im2 (L) is at least a factor of 6 smaller 
than that occurring in z e g(L). 

In Fig. [2] we show estimates of z im (L) and z im2 (L) at p = 0.8. Fits to z + cL^ 2 
give results in perfect agreement with those obtained by fitting t x (L) and z e g(L), see 
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Figure 7. Effective exponents z c g(L) (above) and Zi m (L) (below) versus L w , 
u = 0.29, for the RSIM at p = 0.65, 0.8, 0.85. The dotted lines correspond to the 
final result z = 2.355(16). 

Fig. [31 Note that the data for z im (L) and Zi m 2(L) are very close and provide almost 
equal results. This can be easily explained by noting that 

Zim(L) - z im2 (L) = (r- qU*Jz) bL~"> + • • • , (54) 

where b is a p-dependent coefficient. Since r and q are good approximations of r* and 
q* defined in (|50|) . the prefactor is very small, explaining why the two quantities behave 
identically. It is interesting to note that Z[ m (L) has L~^ 2 corrections which are larger 
than those occurring in z e d(L), see Fig. [2j improved quantities have smaller leading 
scaling corrections but larger subleading ones. 

In the following we only report results for z im (L). Fits of z im (L) at p = 0.8 give 
z = 2.355(8) [z = 2.356(7)] for L min = 16 [L min = 12], where the errors also take into 
account the uncertainty on C/j = 1.500(1). These results vary approximately by ±0.008 
when changing u 2 within [0.74, 0.90], corresponding to the uncertainty on u 2 . We finally 
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obtain the estimate 

z = 2.355(8)[8]. (55) 

Because of the bound f[5"3"j) . the error due to the residual L~ w scaling corrections is 
negligible. This result confirms the one given in (I34|) . 

In Fig. [7] we report z im (L) and z e g(L) for the three values of p we have considered. 
In all cases, the improved exponents are quite close to the final estimate fl55l) . For 
p = 0.85, while z e «(L) was close to 2.25, z- lul {L) is fully consistent with 2.355. As for 
p = 0.65, the difference between z- im (L) and 2.355 is three times smaller than that 
between z e ^{L) and 2.355. The still existing discrepancies can be explained by the 
next-to-leading 0(L~ 2uJ ) scaling corrections, as shown by Fig. [8] where the difference 
Az im (p; L) = z im (p; L) — z irn (p = 0.8; L) is plotted versus L~ 2w . Clearly, Az im (p; L) is 
consistent with zero, for both p = 0.85 and 0.65, if we only consider data with L > 16, 
supporting universality. 

In conclusion, the results of the RSIM provide an accurate estimate of the dynamic 
exponent z, i.e. z = 2.355(16), and a robust evidence of universality, i.e. independence 
on p. 

Finally, we note that the leading amplitude c defined in ( 1261) significantly increases 
with decreasing p. Indeed, we find c ~ 0.03, 0.05, 0.3 for x = 1 and p = 0.85, 0.8, 0.65, 
respectively. 

4-3. Universality of z in the RBIM and ±J Ising model 

We now check the universality of the dynamic exponent z in other RDIs systems, such 
as the RBIM and the ± J Ising model along the paramagnetic-ferromagnetic transition 
line. We first focus on the approximately improved models, the RBIM at p = 0.55 
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Figure 9. The effective exponents z e s(L) and z- lm (L) for the RBIM (rb) at p = 0.55 
and for the ± J Ising model (frb) at p — 0.883, versus L~^ 2 with u>2 = 0.82. 

(the RBIM is improved for p = p* = 0.54(2)) and the ±J Ising model at p = 0.883 
(improvement occurs for p = p* = 0.883(3)). We perform an analysis analogous to that 
presented for the RSIM at p — 0.8, verifying that scaling corrections decay as L L ° 2 , as 
expected on the basis of field theory, and computing for each of them an estimate of z. 
Then, we verify that the results for the other values of p are consistent with universality, 
i.e. that the large observed deviations can be interpreted as scaling corrections. 

In Fig. [9] we plot z eS (L), defined in ([2"T]) . and z im (L), defined in (1521) . versus L~^ 2 
with 0J2 = 0.82. In the case of the ± J Ising model, both z e s(L) and Zi m (L) clearly show 
the expected L~ W2 behaviour. In the case of the RBIM, z- 1Ta {L) shows a clear linear 
trend, while z e s(L) becomes essentially flat as L increases and is close to the RSIM 
estimate z = 2.355(16): indeed, z eS (L) = 2.341(1), 2.336(2), 2.338(4) for L = 16,24,32. 

To verify that the RBIM at p — 0.55 and the ± J model at p = 0.883 have the same 
dynamical critical behaviour as the RSIM at p = 0.8, in Fig. \IU\ we plot the ratio 

r x =i(L) /t x=1 (L) KS iu,p=o.s- (56) 

As L — > 00 the data clearly approach a constant, indicating that all autocorrelation times 
diverge with the same z. The data shown in the figure are well fitted to 6 + 62-^ ~ W2 , with 
b pa 1.74 and b ~ 0.84 respectively for the RBIM and ±J Ising model. 

The results of the fits of r x= i(L), z e s(L), and z- im (L) for the ±J Ising model at 
p = 0.883 are shown in Fig. [TlJ In particular, by fitting t x= i(L) to (|30|) with e = 0.82, 
we obtain z = 2.345(4) and c = 0.0466(8) for L min = 14. The fit of z eS {L) to z + eL^ 2 
gives z = 2.344(6) for L min = 14 and z = 2.342(13) for L min = 20. These results suggest 
the estimate z = 2.345(4) [3]. The error in brackets gives the variation of the estimate 
as u>2 varies by one error bar. By fitting z- im (L) to z + eL~ W2 with 002 = 0.82, we obtain 
z = 2.345(8) (for L min = 14), where the error includes the uncertainty on U%. The 
uncertainty on U2 changes the estimate by ±0.006. 
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Figure 11. Estimates of the dynamic exponent z from fits of T x —i(L) to (|30|) with 
e = 0.82, and of the corresponding z g(L) and Zi m (L) to z + cL~ 82 . The data refer 
to the ± J Ising model at p = 0.883. Some data are slightly shifted along the x axis to 
make them visible. The dotted lines correspond to the final estimate z = 2.345(17), 
see ([SB]). 



As in the case of the RSIM, p* is only known approximately and thus some residual 
leading scaling corrections may still be present. To determine their relevance, we again 
exploit the fact that ratios of amplitudes of leading scaling corrections are universal, 
and the bound 

\co,n(p = 0.883)/co,u(p = 0.9)| < 1/5. (57) 

The error on the estimate obtained from z e s(L) due to possible residual scaling 
corrections can be estimated as in the case of the RSIM at p — 0.8, obtaining ±0.015. 
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This error is significantly smaller when z lWL {L) is considered: the estimate varies by 
±0.003. In conclusion the most precise estimate of z for the ±J Ising model is obtained 
by using z- im (L). We quote 

z = 2.345(8)[6]{3}, (58) 

which is in good agreement with the RSIM result (155]) . 

Let us now consider the RBIM at p — 0.55. By fitting z- 1Ui (L) we obtain z = 2.336(7) 
for L min = 12, and z = 2.335(9) for L min = 16. We obtain z = 2.335(9) [4], where the 
error in bracket gives the change in the estimate as u>2 varies by one error bar. As in the 
case of the RSIM and the ±J Ising model, we must also estimate the error due to the 
residual scaling corrections. Using the results reported in [6], we find that these 
corrections can at most change the estimate of z by ±0.005. Our final result is therefore 

z = 2.335(9)[4]{5} . (59) 

In Fig. [T2l we show z e g(L) for other values of p, i.e. for the RBIM at p — 0.7 and 
for the ±J Ising model at p = 0.83,0.90. They are plotted versus L~ u , which is the 
expected leading scaling corrections. As it was observed for the RSIM, see Fig. [TJ the 
results appear strongly p-dependent and it is not clear from the data that z e s(L) has a 
model- and p-independent limit as L — > oo. In any case, we can show that these data 
are still consistent with universality if the expected scaling corrections are taken into 
account. Let us again consider the difference 

Az RS (p; L) = z eS (p- L) - z eff (RSIM,p = 0.8; L), (60) 

which, as discussed in the preceding section, should behave as 

Az RS (p; L) « e!(p)L- w + e 2 L" e , (61) 

with e = 2u,U2 for L — > oo. We recall that e"i ~ en, cf. ( 1281 . In order to show 
consistency with universality, we fit Az R s(p; L) to the ansatz ( 1611 . determining en. 
Then, we consider Z7 2 2 and determine the leading correction-to-scaling amplitude c 2 2,n 
defined in (|4"2"]1 . Finally, we verify that the ratio S22 = en/c22,n is independent of p and 
of the model, and that it agrees with the RSIM estimate (jlB"]) . s 2 2 = 9(2). 

As shown in Fig. [13], good fits of Az R s(p; L) to ( TBTj) are obtained by taking e = 2u. 
They give e"i = -1.0(2) for the RBIM at p = 0.7, and e x = 1.2(2), -0.5(1) for the 
± J Ising model at p = 0.83, 0.9, respectively. The amplitude 022,11 can be estimated 
analogously. We obtain c 22 ,ii = -0.17(3), 0.10(2), -0.05(1), respectively for the RBIM 
at p — 0.7, and the ± J Ising model at p = 0.83, 0.9. These results give 

s 22 = 6(2), 12(3), 10(3), (62) 

which are in substantial agreement with the estimate f ]46l) obtained from the RSIM. 
These results fully support the interpretation of the anomalous behaviour of the data 
shown in Figs. [12] and d3] as an effect of scaling corrections. 

In Fig. [12] we also show z im (L). For the ±J model at p — 0.83, the improved 
estimator is significantly closer to z ~ 2.35 than z e s{L). In the two other cases 
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Figure 12. Effective exponents z e g(L) and z- lm (L) for the RBfM at p — 0.7 (rb) and 
the ±J Ising model at p = 0.83,0.90 (frb). The dotted lines correspond to the final 
estimate z — 2.35(2). 
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deviations are still large, though the data seem to approach faster the limiting value 
z ~ 2.35 obtained by considering the improved models. In Fig. [THwe show the difference 
z im (L) — 2.35 versus L~ 2uj , which would be the leading scaling correction if z im (L) 
were exactly improved. The data for the ± J Ising model converge to zero, confirming 
universality On the other hand, the RBIM data apparently extrapolate to a slightly 
positive value. However, if we include an additional correction term (either an L~ u term, 
since improvement is only approximate, or an L~ W2 term), the data are again perfectly 
consistent with universality 
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Figure 14. The difference z im (p; L) - 2.35 versus i^ 2tJ with w = 0.29 for the RBIM 
(rb) and the ± J Ising model (frb). 



4-4- Summary of the results for the equilibrium relaxational dynamics atT c 

In this section we have studied the Metropolis equilibrium dynamics in the RSIM, the 
RBIM, and the ±J Ising model along the paramagnetic-ferromagnetic transition line, 
which belong to the same static universality class. We have verified that the exponent z 
is the same for the RSIM, the RBIM, and the ± J Ising model for values of the disorder 
parameter p that make these models approximately improved. We have obtained 
z = 2.355(16), z = 2.335(18), and z = 2.345(17) respectively for the RSIM at p = 0.8, 
the RBIM at p = 0.55, and the ± J Ising model at p — 0.883. For the other values of p we 
have not been able to determine z with comparable precision. We have however verified 
that the dynamic behaviour is always consistent with universality once the expected 
scaling corrections are taken into account. In the analyses we have presented, scaling 
corrections play a very important role. We have explicitly verified the FT prediction 
that dynamics does not introduce new RG irrelevant operators that are more relevant 
than the two leading ones occurring in the statics. Therefore, scaling corrections are 
characterized by the same universal exponents that control the nonasymptotic behaviour 
in static quantities, i.e. u = 0.29(2) and u 2 = 0.82(8). As a consequence, the leading 
L~ kuJ scaling corrections are absent in dynamic observables at the same value p* of the 
disorder parameter p determined by considering static quantities. 

Once universality has been checked, we can use our results for the RSIM, the RBIM, 
and the ±J Ising model, to obtain a final estimate for z. We quote 

z = 2.35(2), (63) 

which includes all results obtained in the previous sections. 

Our result ( 1631) significantly improves earlier MC estimates [HI [19J, [23] of z obtained 
in equilibrium MC simulations. Reference [19J considered several values of p in the range 
1 > p > 0.6. The final estimate z = 2.4(1) was essentially derived from the data at 
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p = 0.8, where the finite-size behaviour appeared least dependent on L. The different 
values of z obtained for other values of p were interpreted as a crossover phenomenon. 
A smaller estimate z = 2.2(1) was found in [23], by MC simulations at p = 0.85. This 
may be explained by the effect of scaling corrections, which, as shown by Fig. [TJ give 
rise to an apparent smaller value of z if they are not taken into account. 



5. Off-equilibrium estimate of the dynamic critical exponent z 

The exponent z can also be determined by performing off-equilibrium simulations 
[30~t I3TI [T6] . One starts from a disordered configuration and observes the relaxation at 
T c on sufficiently large lattices. In the infinite-volume limit the magnetic susceptibility 
X is expected to increase with the MC time t as 

x(t) = c t" (i + c n r ul + c l2 r 2 ^ + ■■■ + c 21 r v > + •••), (64) 

where 

2 — 77 

P = (65) 

Using the estimate obtained in Sec. 15 z — 2.35(2), and [6J r\ = 0.036(1), we predict 
p = 0.836(7). Moreover, according to the FT perturbative analysis [301 ED fTB"] . 
the leading scaling-correction exponents should be the same as those that occur in 
equilibrium (static or dynamic) correlation functions. Therefore, we expect 

v x = - = 0.123(9), v 2 = — = 0.35(3), (66) 

z> z 

where we have used u = 0.29(2) and u 2 = 0.82(8). Moreover, the leading scaling 
correction proportional to t~ Vl (and also all corrections of the form t~ kvi ) vanishes in 
improved models. 

Equation fl64|) is valid only in the infinite- volume limit. For a finite system of size 
L we expect 

X (t,L) = C t" E (tL- z )(l + C lx t^E x (tL- z ) + • ■ ■) (67) 

where Ei(x) are universal functions satisfying Ei(0) = 1 and Eq(x) ~ x~ p , E x (x) ~ x vi 
for large x. 

The off-equilibrium relaxational dynamics of the RSIM has already been 
investigated in [20, 22J for various values of p in the range 1 > p > 0.4. Their results do 
not agree with the above- reported predictions. References [201 122] obtain z = 2.62(7) 
and z = 2.6(1), respectively, independently of the dilution parameter p. They also 
estimate the leading correction-to-scaling exponent v x . Given their estimate of z, this 
allows them to estimate to. They quote uj = 0.50(13) and uj = 0.61(6), respectively. It 
is quite difficult to reconcile these results with the FT predictions; in particular, the 
absence of corrections proportional to t~ 012 , i.e. controlled by the leading exponent 
to = 0.29(2), is quite surprising. Another numerical MC work [21] investigated the 
nonequilibrium relaxation dynamics of the ± J Ising model, and showed an apparent non- 
universal dynamical critical behaviour along the paramagnetic-ferromagnetic transition 




Figure 15. Effective exponent p e g(t,L) for L = 64,96, 128. Here t is the number of 
MC sweeps. 



line. Also these results are in contrast with the FT predictions reported at the beginning 
of section. 

In the following we further investigate this issue. We study the Metropolis dynamics 
of the RSIM at p = 0.8 after a quench from T = oo to T c . This represents a nontrivial 
check of the FT predictions, since the estimates (165]) and fl66|) are quite precise. Since 
the model is approximately improved [p* = 0.800(5)], we expect that C±k vanishes for 
all values of k, and thus we predict 

= C t p (1 + C 21 t-" 2 + ■••), (68) 
p = 0.836(7), v 2 = 0.35(3). 

As in the equilibrium case, we define an effective exponent 

w(t) . ^Xi, (69) 
which behaves as 

Pcfr(t) =p + ct~ V2 + ... (70) 

for t —>■ oo. On a finite lattice, flTOj) is replaced by 

PeS (t,L) =p + e (tL- z ) (71) 

where we have neglected large-t (scaling) corrections and Cq{x) is a universal function 
(apart from a normalization of the argument) such that eo(0) = and Cq{x) — > —p for 
x — > oo. 

We have performed off-equilibrium MC simulations on lattices of size L = 64, 96, 128 
at (3 = 0.2857430 [our presently best estimate of (3 C is (3 = 0.2857431(3)]. For each 
lattice size we average over N s = 320000 disorder configurations. For each disorder 
configuration we start from a (different) T = oo configuration and perform 2000 
Metropolis sweeps, using the algorithm described in Sec.l4.1land|Appendix B| In Fig. [T5l 
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Figure 16. The effective exponent p e ff(t,L) versus t Tesc = t(L/128) 



we show p c g(t,L) for L = 64,96,128. It clearly shows finite-size corrections, and, for 
each L, p e ff(t,L) follows the infinite-volume curve up to an L-dependent value t mgcK (L). 
As shown by Fig. CEl where p e g(t,L) is plotted versus t resc = t(L/128)~ z , finite-size 
effects are consistent with flTTj) . Thus, the value t max (L), after which finite-size effects 
cannot be neglected, increases as L z . Infinite- volume quantities, such as p e s(t), must 
be obtained from the data at t < t max (L). Fig. [T5l indicates that, with the statistical 
errors of our data, t max (L) ps 120,600 for L ps 64, 128. Since p e s(t,L) is defined using 
data at t and 2t, this implies that, for L = 128, only data corresponding to t < 1200 
have negligible finite-size effects within our error bars. Finite-size effects give rise to a 
systematic error in the estimate of p. As is clear from Fig. [T5], they yield smaller values 
of p, and therefore larger values of z. 

In Fig. [17] we plot p e g(t,L) for L = 64, 128 versus t~ V2 with v 2 = 0.35. Finite-size 
effects are negligible for t~ V2 > t m3iX (L)~ V2 m 0.17,0.10, for L = 64,128, respectively. 
The data satisfying this inequality clearly follow a unique curve, which is expected to 
behave as p + ct~ V2 for sufficiently large values of t. The data plotted in Fig. [17] clearly 
show such a behaviour in the region t~ V2 < 0.4 (corresponding to t > 10), and are 
perfectly compatible with the values p = 0.836(7) and v 2 = 0.35(2). This is also shown 
by Fig. [TBI where we plot the results of fits of p e s(t, L) for t min < t < t max to p + dr V2 
with v 2 = 0.35, for L = 64, 128. 

The above results provide a nice check of the results of the previous section 
and confirm the RG prediction that the off-equilibrium relaxational critical dynamics 
is characterized by the same dynamic exponent z and the same scaling-correction 
exponents uj and uj 2 as the equilibrium critical dynamics. Note that the results of 
Fig. [171 rule out the larger estimates of z obtained in [201 122], z = 2.62(7) and z = 2.6(1), 
corresponding to p = 0.750(20) and p = 0.755(29), respectively (using [6] rj = 0.036(1)). 
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Figure 17. Effective exponent p c s(t, L) versus t V2 with v 2 = 0.35. Finite-size effects 
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Figure 18. Results of fits of p e s(t, L) for t mm < t < t max to p + ct~ V2 with v 2 — 0.35, 
for L = 64, 128. The dotted lines correspond to the prediction p — 0.836(7) obtained 
by using the equilibrium result z = 2.35(2) and 77 = 0.036(1) [B]. 

6. Conclusions 

In this paper we have studied the purely relaxational dynamics (model A) in randomly 
diluted Ising (RDIs) systems. According to standard RG arguments applied to 
dynamics, the dynamic critical behaviour in such systems should belong to the same 
model-A dynamic universality class. If this description is correct, the dynamic exponent 
z is the same in all RDIs systems and the leading scaling corrections are controlled 
by the same RG operators that appear in the statics and therefore are characterized 
by the static correction-to-scaling exponents to = 0.29(2) and 002 = 0.82(8). For the 
same reasons, in the case of improved Hamiltonians, leading scaling corrections should 
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also be absent in dynamical quantities. Therefore, improved models are expected to 
provide the most precise estimates of universal dynamic quantities. For instance, in 
FSS studies at the critical point, corrections to scaling decay as L~°- 82 in improved 
models, while in generic RDIs systems the approach to the infinite-volume limit is much 
slower, corrections decaying as L -0 29 . 

The main results of our analysis can be summarized as follows. 

(i) We consider the RSIM at p = 0.8, the RBIM at p = 0.55, and the ± J Ising model 
at p = 0.833, at the critical point. These three models are approximately improved 
(the best estimates of p*, the value of the disorder parameter corresponding to an 
improved model, are 0.800(5), 0.54(2), 0.833(3) in the three models, respectively 
[61 [7j). We perform high-statistics equilibrium MC simulations on lattices L 3 , 
L < 64, using the Metropolis algorithm (for the RSIM and the RBIM a small 
modification is needed to ensure ergodicity, see Appendix B ). We determine the 
exponent z, obtaining z = 2.355(16) for the RSIM, z = 2.335(18) for the RBIM, 
and z = 2.345(17) for the ± J Ising model. These results are in perfect agreement, 
providing strong support to the FT prediction that all RDIs models belong to the 
same dynamic model-A universality class. We also investigate in detail the scaling 
corrections: they are perfectly consistent with a behaviour of the form L~ W2 , with 
u>2 = 0.82(8). Again this is in agreement with the FT analysis. Our final result is 

z = 2.35(2). (72) 

Note that, while we confirm the general scenario predicted by field theory, there is 
a quantitative difference between our result and that obtained by resumming the 
perturbative expansions at two and three loops, z ~ 2.18 [281 [29]. This may be due 
to a poor convergence of the perturbative FT expansions. The apparent agreement 
with the 0(y/e) result [24J z = 2 + ^/6e/53, which would give z ~ 2.336 for e = 1, 
is likely only by chance. 

(ii) We investigate the Metropolis dynamics in equilibrium in the RSIM, the RBIM, 
and in the ± J Ising model for other values of p. Here, as expected, corrections are 
very strong. In the FSS analysis, the leading term is expected to decay as L~ w , 
us = 0.29(2). We are not able to determine z in these models as accurately as 
in improved models. In any case we verify that the L-behaviour of the MC data 
at T c is consistent with universality and with the constraints imposed by the RG 
(universality of ratios of correction-to-scaling amplitudes). 

(iii) We have no evidence of two different universality classes depending on the disorder 
strength [441 145]. In particular, we show that the critical behavior is not influenced 
by the geometrical structure of the vacancies and does not depend whether the 
vacancies percolate or not. Indeed, since site and bond vacancies percolate for 
p < 0.688 and p < 0.751, respectively, in the improved RSIM (p = 0.8) vacancies 
form finite clusters, while in the improved RBIM (p = 0.55) vacancies percolate. 
Nonetheless, the critical behavior is the same. 
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(iv) The results for the ±J Ising model show that frustration is irrelevant for the 
relaxational behaviour along the paramagnetic- ferromagnetic transition line. It 
is worth mentioning that this is not true for the cluster dynamics. In that case the 
exponent z in the ±J Ising model is much larger than in the RSIM and RBIM. 
In the frustrated case we obtained z ~ 1.6 [7], while in the second one simulations 
indicate z < 0.5 [23] . 

(v) We investigate the off-equilibrium relaxational dynamics in the RSIM at p = 0.8. 
We start from disordered T = oo configurations and observed the relaxation 
at T = T c . The results show that our equilibrium estimate z = 2.35(2) is 
perfectly consistent with the off-equilibrium MC data. In the analysis particular 
care has been taken to avoid finite-size corrections. Moreover, the large-time 
scaling corrections are consistent with what is expected on the basis of field theory 
[301 [31j [16] . Since the model is improved, we do not observe corrections proportional 
to t~ w l z \ instead our data show corrections that are proportional to t~ w ' 2 l z . Here cu 
and oji are the static correction-to-scaling exponents, lj = 0.29(2) and u 2 = 0.82(8). 

The relaxational critical behaviour within the RDIs universality class is also relevant 
for the so-called model-C dynamics, where the order parameter couples with a conserved 
secondary density [llj. In the case of the 3D RDIs universality class, the asymptotic 
critical behaviours of the model-A and model-C dynamics are described by the same 
stable fixed point. Therefore, they are expected to have the same dynamic exponent z. 
This is essentially related to the fact that the specific-heat exponent of RDIs systems, 
a = —0.049(6), is negative [46] . A drastic change occurs in the approach to the 
asymptotic behaviour, which is expected to be much slower in model C |47j. The 
coupling with a conserved scalar density gives rise to very slowly decaying 0(£ _Wc , L _tJc ) 
scaling corrections with0 

3 1a.. , 
= 2 — v = ~Yu = °- 036(4) ' (73) 
which is much smaller than the leading scaling-correction exponent of the model-A 
dynamics, which is to = 0.29(2). 
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(j When the specific-heat exponent a is negative, the asymptotic critical behaviour of model C is 
the same as that of model A, because they have the same stable fixed point [46], [11]. Nevertheless, 
the coupling between the order parameter ip(x) and the conserved scalar density e(x) introduces a 
new irrelevant RG perturbation, which is not present in the model A and which gives rise to very 
slowly decaying scaling corrections. The RG dimension y c of the Hamiltonian coupling term Ti ve = 
7o / d d xe(p 2 can be computed by using nonperturbative scaling arguments: y c = y ip 2+y £ —d = l/v—d/2. 
This implies that there are 0(£~" c ) scaling corrections to the asymptotic critical behaviour, with 
uj c = —y c . Using [6] v = 0.683(2), one obtains uj c = 0.036(4). 
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Figure Al. Estimates of the leading correction-to-scaling exponent cj. The dotted 
lines correspond to our final result cu = 0.29(2). 



Appendix A. Estimate of the leading correction-to-scaling exponent lu 

In this appendix we compute the leading correction-to-scaling exponent lu. We use 
the method discussed in j6] and combine the data of [6J with those obtained here. 
We consider the quartic cumulants f7 22 and U d , cf. (14TI) . at fixed £/L = 0.5943, at 
p = 0.85, 0.8.0.65 on lattices of size L < 192. 

As in [6], in order to estimate lu, we consider the differences 

A 22 (pi,P 2 ; L) = U 22 ( Pl ; L) - U 22 (p 2 ; L), (A.l) 
A d ( Pl ,p 2] L) = U d ( Pl ; L) - U d (p 2 ; L). (A.2) 

Universality implies that 

A ss c AA1 L-" + c A ,i2L- 2u + ■■■ + c A , 2 i^^ 2 + • • • (A.3) 

In the case of A 22 , fits to cL~ w provide good and stable results. In the case of A^, an 
additional correction term is needed in order to obtain an acceptable y 2 . Hence, we fit 
A d to C\L~ W + c 2 L~ £ . In Fig. lAll we show the results as a function of L m i n , the minimum 
lattice size used in the fit. They become independent of p%, p 2 , and of the considered 
quantity as L min increases. The most stable results are obtained by taking pi = 0.85 
and p 2 = 0.65. From the results shown in Fig. IA.ll we obtain the estimate 

u = 0.29(2). (A.4) 

This estimate is more precise than previous ones and is consistent with lu = 0.33(3) 
obtained in [6] by using part of the data at p — 0.65 and p = 0.8. For comparison, the 
FT six-loop result [9] is lu — 0.25(10) (we also mention the five-loop result lu = 0.32(6) 

of ma). 
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Appendix B. Metropolis algorithm for randomly diluted Ising models 

We have implemented the standard local Metropolis algorithm with the acceptance rate 

A = mm[l, exp(- (3[H' -H])], (B.l) 

where "H' and TC correspond to the Hamiltonian evaluated for the proposal and for the 
given spin configuration, respectively. The proposal is generated by flipping the sign of 
the spin at a single site x of the lattice. Hence Ji! — 7i depends only on the values of 
the spins at the site x and at its neighbours y. 

To speed up the simulation we use multispin coding (see, e.g., [IB]), evolving in 
parallel nut systems (nbit = 64 in our case). For each of them we use a different set of 
disorder variables. The implementation in the RSIM and RBIM is more complicated 
than in the standard Ising model and the ± J Ising model, since the sum over the nearest 
neighbours assumes one of the 13 values {—6, —5, 5, 6} and not only the 7 even values 
{-6, -4,. ..,4,6}. 

We perform the single- site update sequentially, moving from one site to one of its 
neighbours in a typewriter fashion. This causes problems with ergodicity. This can be 
understood by considering the isolated lattice sites, i.e. the sites x such that p y = 
(RSIM) or J <xy> = (RBIM) for all neighbours y. For an isolated site the difference 
H! — 7i always vanishes, so that, using the acceptance rate (IB. 1 j) . the proposal is always 



accepted. Hence, going through the lattice twice, the spins on the isolated sites go 
back to their values. Therefore, the configuration restricted to the isolated sites is not 
changed. Note that the problem is not restricted to isolated sites only. For example, for 
the one-dimensional chain one can easily prove that a regular update sweep using the 
acceptance rate fIB.lj) is not ergodic. There are many ways to avoid this problem. For 



performance reason, we prefer to update the spins sequentially. To avoid the problem 
the spin flip is proposed with probability w strictly smaller than one. We have chosen 
w = 0.9. Note that the problem occurs only in the RSIM and in the RBIM. For the ±J 
Ising model the standard Metropolis update can be used. 
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